#ifndef __DFT_HARD_SPHERE_H
#define __DFT_HARD_SPHERE_H

#include<dftFFT.h>

#define NUMBER_TO_MOLAR 1.6606

#define DFT_DEBUG_HARD_SPHERE_n 3
#define DFT_DEBUG_HARD_SPHERE_WHAT2 3
#define DFT_DEBUG_HARD_SPHERE_WHAT2_DETAIL 4


// Note: '2' as in 'rho2' or 'rho2Hat' or 'WHat2' refers to vectors living on reflected DAs 
// (twice the size in each dimension of the "physical" DAs.
// Likewise for the actual DAs.

typedef struct DFTHardSphere : public DFTFFT {
  PetscTruth     check_pf;
  PetscInt       numAlpha;
  PetscInt       dims2[3];
  DA             rhoDA2;
  DA             singleDA, singleDA2;
  DA             weightsDA, weightsDA2;
  Mat            F2;       // Discrete Fourier Transform operator for the DA
  /* Work Vecs */
  /** Rho DAs **/
  Vec           rho2, mu2;
  /** Weights DAs **/
  Vec           n2, PhiDer2, PhiDer2Hat;
  /** Single DAs **/
  Vec           z2; // A purely auxiliary accumulator Vec.
  /* Work Vec sets: split by species or by alpha */
  Vec           *rho2BySpecies,  *rho2HatBySpecies;
  Vec           *mu2BySpecies,   *mu2HatBySpecies;
  Vec           *WHat2BySpecies, *WHat2SpeciesByAlpha;
  Vec           *n2ByAlpha,      *n2HatByAlpha;
  Vec           *PhiDer2ByAlpha, *PhiDer2HatByAlpha;
} DFTHardSphere;

typedef struct {
  PetscScalar w[6];
} Weights;




EXTERN PetscErrorCode DFTHardSphereCreate(DFTHardSphere *dftHS);
EXTERN PetscErrorCode DFTHardSphereDestroy(DFTHardSphere *dftHS);
EXTERN PetscErrorCode DFTHardSphereCalculateMu(DFTHardSphere *dftHS, Vec rho, Vec muHS);
//EXTERN PetscErrorCode CalculatePressure(const PetscInt numSpecies, const PetscReal R[], const PetscTruth isConfined[], const PetscReal rhoBath[], PetscReal *P);
//EXTERN PetscErrorCode CalculateRhoBath(DM dm, const PetscInt numSpecies, Vec rhoBathVec);
//EXTERN PetscErrorCode DFTCalculateMuHSBath(const PetscInt numSpecies, const PetscReal R[], const PetscTruth isConfined[], const PetscReal rhoBath[], PetscReal muHSBath[]);

#endif /* __DFT_HARD_SPHERE_H */
